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Abstract 

Several research groups have recently reported ab initio calculations of the melting properties 
of metals based on density functional theory, but there have been unexpectedly large disagree- 
ments between results obtained by different approaches. We analyze the relations between the 
two main approaches, based on calculation of the free energies of solid and liquid and on direct 
simulation of the two coexisting phases. Although both approaches rely on the use of classical 
reference systems consisting of parameterized empirical interaction models, we point out that 
in the free energy approach the final results are independent of the reference system, whereas 
in the current form of the coexistence approach they depend on it. We present a scheme for 
correcting the predictions of the coexistence approach for differences between the reference and 
ab initio systems. To illustrate the practical operation of the scheme, we present calculations 
of the high-pressure melting properties of iron using the corrected coexistence approach, which 
agree closely with earlier results from the free energy approach. A quantitative assessment is 
also given of finite-size errors, which we show can be reduced to a negligible size. 

1 Introduction 

The study of solid-liquid equilibrium by computer simulation has a long history, going back to 
the classic work of Alder and Wainwright [[ij on the hard-sphere system in the 1950's. Many 
techniques have been used to determine the pressure-temperature relation at equilibrium and 
other melting properties, such as the volume and entropy of fusion. In the last few years, there 
has been an upsurge of interest in the accurate ab initio treatment of the melting properties 
of real materials [||, |3|, ||, [5|, §, 0, g, which has focused attention again on the techniques 
used to locate the melting transition. Recent ab initio work has been based on two main 
approaches. The first locates the melting transition by requiring equality of the Gibbs free 
energies, which are calculated ab initio for liquid and solid || [||, |], ||, [| ; we call this the 'free 
energy' approach. The second proceeds by fitting a potential model to ab initio calculations 
and using this to simulate a system containing liquid and solid in coexistence 0, |8| |9) ; we call 
this the 'coexistence' approach. If appropriate measures are taken, the two approaches should 
clearly give the same results. The purpose of this paper is to analyze the relation between the 
two approaches and to propose what these 'appropriate measures' should be. We illustrate our 
analysis by presenting ab initio calculations on the high-pressure melting of iron performed 
using the coexistence approach, which we compare with earlier free-energy results from the 



same ab initio technique 0]. We shall show that, once the necessary corrections have been 
applied, the two approaches yield the same results. 

The free-energy approach has been well established for many years for calculations based 



on classical interaction models (see e.g. Refs. [10, 11, [L2|]). Typically the procedure for the 



solid has been to start from the harmonic approximation at low temperatures; the free energy 
of the high-temperature anharmonic system is then obtained by using the Gibbs-Helmholtz 
relation to integrate the internal energy given by molecular dynamics (MD) simulation [10, |ll] [. 
Alternatively, the anharmonic free energy has sometimes been obtained by thermodynamic 
integration starting from a reference model such as the Einstein solid For the liquid, 

a common procedure has been to obtain the free energy at one thermodynamic state from 
the work done in reversible expansion [10, [l^] to low density or heating to high temperature. 
The Gibbs-Helmholtz relation is then used to obtain the free energy at other states. We also 
mention an important alternative approach, known as 'Gibbs-Duhem' integration, which allows 
the boundary between coexisting phases to be directly mapped out |[^, Q. 

In early work, the interaction models were parameterized by fitting to experimental data, 
but recent years have seen a major shift towards the calculation of melting properties |2|, ||, 
13 H' @' S @ f rom a b initio methods based on density- functional theory (DFT) | |15| |. In DFT, 
the total energy function of a system is determined by the approximation used for exchange- 
correlation energy E xc . An important ambition then follows: the determination of melting 
properties with no statistical-mechanical or other approximations except those inherent in E KC 
itself. This raises major new issues, because it is extremely costly to perform ab initio MD 
(AIMD) simulations on large systems long enough to reduce statistical errors to an acceptable 
level. The cost is particularly great for crystals, since extensive electronic Appoint sampling may 
be needed. Long accepted methods like reversible expansion to low density are impracticable 
with AIMD. Instead, it is essential to use empirical interaction models that closely mimic the 
DFT system. The usual techniques are then employed to treat these models, and in some ab 
initio work thermodynamic integration is used to obtain the free energy difference between the 
model and DFT systems [Q, ||, |], ||, [|. In parallel with these developments, the advantages of 
avoiding intricate free-energy calculations by using the coexistence approach have made this 
route popular ||, M. Again, the simulations must be based on empirical models matched to 
DFT data. 

The present work was stimulated by recent reports of ab initio or nearly ab initio cal- 
culations on the melting of a number of materials, including Si [§], Al ||, [f|, Cu [16] and 
Fe j|, [|, ^, [8| . Our own work on Al ||, || and Fe Q ||| was based on the free-energy approach, 
which used an inverse-power reference model, with thermodynamic integration to calculate the 
difference between the ab initio and reference free energies. A major effort was made to ensure 
that finite-size and other systematic errors were reduced to an almost negligible level. The 
other ab initio work on Fe, by two independent groups, used the coexistence approach [[?], ||]. 
The results for the melting curve differed substantially from each other and from our results, 
and we urgently need to understand the reasons for the disagreements. We aim to shed light 
on possible reasons here. 

We shall explore a number of technical issues. The first concerns the correction of the 
coexistence method for errors due to the difference between the ab initio system and the 
empirical model that mimics it. A second issue is the fitting of models to ab initio data. We 
have already discussed this in depth for the free-energy approach, so here we focus mainly on 
coexistence. We want to study what physical quantities should be fitted, and how to tell if 
the fit is good enough. A third important issue concerns finite-size errors, which arise mainly 
from the limited system sizes that can be handled ab initio. It has been claimed |16] that 
the coexistence approach suffers less from size errors than the free-energy approach. We shall 
demonstrate that in fact this type of error affects both approaches in essentially the same way. 
One thing we shall not do is to pass judgment on which method is 'better', since we shall argue 
that the answer depends on what one wishes to achieve, and that both approaches are vital. 



However, we shall comment on the strengths and weaknesses of both. 

In the next Section, we define the problem, summarize briefly the free-energy and coex- 
istence approaches as they have been applied in practical ab initio calculations, analyze the 
corrections that need to be applied in the two approaches, and discuss size errors. Section |j| 
then presents our new calculations on the high-pressure melting properties of Fe using the coex- 
istence approach. We shall compare with our earlier free-energy results and show the practical 
necessity of the corrections outlined in Sec. ^. The final Section gives further discussion and 
summarizes our conclusions. 

2 Theory of the two approaches 

2.1 Definition of the problem 

We start by defining the problem to be addressed. At a given pressure p, two phases are in 
thermodynamic equilibrium when their Gibbs free energies G(p, T) are equal. We regard G as 
derived by the relation G = F + pV from the Helmholtz free energy F(V,T). The ab initio 
value of F at any volume and temperature, denoted by F\i(V, T), is given in classical statistical 
mechanics by: 

F AI (y,T) = -fe B Tln|^i^|dr 1 ...dr we -^Ai(r 1 ,..r J v)| ; (1) 

where £/Ai( r i> ■ ■ ■ *~n) is the ab initio total energy as a function of the positions ri, . . .r/y of the 
iV nuclei, which range over the volume V of the cell containing the system. For the purpose 
of treating phase equilibria, we need the free energy per atom Fai/N in the thermodynamic 
limit N — > oo, V — > oo at constant number density N/V . 

Even though the coexistence approach does not work directly with free energies, we shall 
show later that the errors that need to be overcome in this approach, as well as in the free- 
energy approach, can be formulated as free-energy errors. But the errors in calculating Faj are 
of two kinds: electronic-structure errors, i.e. imperfections in the calculation of the ab initio 
total energy Uai at each set of nuclear positions; and statistical-mechanical errors, i.e. errors 
in calculating Fai from the given Uai, and in taking the thermodynamic limit. 

For the purposes of this work, electronic-structure errors are irrelevant: our sole concern is 
the treatment of phase equilibrium using some given algorithm for computing Uai- The 'free 
energy' and 'coexistence' routes differ only in the way they address the statistical mechanics. 
Our problem is therefore to assess and compare the ways that the two approaches control the 
purely statistical-mechanical errors. When we come to the practical calculations in Sec. ||| we 
shall compare the results of the two approaches applied using exactly the same algorithm for 
U M - 

2.2 The free-energy approach 

In the free-energy approach [||, f|, ||, [l7| , we use a reference system with total energy function 
f7 re f (ri, . . . rjv), whose solid and liquid free energies F re f are calculated for very large systems, 
so that finite-size errors are negligible. Then the only demanding problem is the calculation of 
the difference AF = Fai — F re { , which is accomplished by thermodynamic integration (l^] : 

AF= f 1 (AU)\ dX , (2) 
Jo 

where AU = Uai — U TC {, and the thermal average (-)\ is taken in the ensemble generated by 
the total-energy function U\ = (1 — A)f7 re f + XUai- In practice, (AU}\ is computed as a time 
average in an AIMD simulation whose dynamics is governed by U\. The main kinds of error 
are: statistical errors in the evaluation of (AU)\; integration errors due to inadequate numbers 



of A points; and finite-size errors. In the systems studied so far by this approach, all these 
errors can be brought under tight control. 

It is a crucial feature of this approach that F^j, calculated as F re f + AF, does not depend 
on the choice of reference system, provided all technical tolerances are set so as to suppress the 
errors we have just mentioned. However, for reasons expounded in detail elsewhere j|, 
the choice of C/ re f has a major influence on the overall computational effort, which is minimized 
by reducing as much as possible the strength of the fluctuations {SAU 2 ) X /N, where 5AU = 
AU — (AU)x- In particular, minimization of this fluctuation strength is important in ensuring 
that -F re f accounts for almost all the free energy, so that size errors in the small residue AF are 
negligible. This can be achieved by using a parameterized reference model whose parameters 
are adjusted so as to minimize (5AU 2 )\/N, as in our work on the melting of Fe. In this 
approach, it is an advantage that different reference systems can be used for solid and liquid, 
since for many materials it may be difficult to create ?7 re f functions that mimic Uai with high 
precision in both phases. 



2.3 The coexistence approach 

The coexistence approach is also based on a reference model that mimics the ab initio system. 
(This model is sometimes called by names such as 'optimised potential model' @ 3 19|, but here 



we give it the same name as in the free-energy approach.) Various ways have been used to fit 
the reference model to the ab initio system. Since it will be relevant later, we note that recent 
coexistence work on Fe has used the 'force-matching' procedure of Ercolessi and Adams in 
which the reference parameters are adjusted so that the ab initio atomic forces are reproduced 
as well as possible for representative sets of atomic positions. 

There are also several ways of using the reference model to simulate coexisting phases, 
and hence to determine the phase boundary between them. In the work of Morris et al. |J, 
coexisting solid and liquid Al were simulated with the total number of atoms N, volume V and 
internal energy E fixed. They showed that, provided V and E are appropriately chosen, the 
two phases coexist stably over long periods of time, and the average pressure p and temperature 
T in the system give a point on the melting curve. The underlying concept is that the mean 
volume per atom v = V/N is given by: 

v = (1 - x l )v s {p) + xiviip) , (3) 

where xi is the fraction of the atoms in the liquid phase, with v s (p) and vi{p) the volumes per 
atom in the coexisting solid and liquid as a function of pressure. For fixed v, the pressure p 
traverses a certain range as x\ goes from to 1. As x\ varies in this way, the mean internal 
energy per atom e = (1 — xi)e s (p) + xiei(p) also traverses some range (here, e s (p) and ei(p) are 
the internal energies per atom in the two phases). Provided e lies in this range for the given 
v, the simulation will yield stably coexisting solid and liquid. An alternative procedure would 
be to simulate at constant (iV, V, T). Then coexistence will be achieved for a given v provided 
T is chosen so that the corresponding p on the melting line curve falls in the range specified 
by Eq. (||). Yet another approach was used in the work of Laio et al. [0] on the high-pressure 
melting of Fe; this used constant-stress simulations, with enthalpy almost exactly conserved. 
The approach of Belonoshko et al. H is different again. Here, the (N,p,T) ensemble is used. 
The system initially contains coexisting solid and liquid, but since p and T generally do not lie 
on the melting line, the system ultimately becomes entirely solid or liquid. This approach does 
not directly yield points on the melting curve, but instead provides upper or lower bounds, so 
that a series of simulations is needed to locate the transition point. Whichever scheme is used, 
some way is needed of monitoring which phases are present. In the (TV, V, E) method of Morris 
et al. ||, graphical inspection of particle positions appears to have been used, supplemented 
by calculating of radial distribution functions to confirm the crystal structure of the solid. In 



the (N, p, T) method of Belonoshko et al. || , the primary diagnostic is the discontinuity of 
volume as the system transforms from solid to liquid. 

The coexistence calculations presented here on the high-pressure melting of Fe employ the 
(N, V, E) method applied to reference systems consisting of the 'embedded-atom model' 
fitted to ab initio data. For geophysical reasons, we are interested in pressures near that at the 



boundary between the Earth's inner and outer cores, namely 330 GPa [21]. In this pressure 



region, the most stable crystal structure just below the melting curve is believed to be hexagonal 



close packed (h.c.p.) p2| , 23], and we assume here that melting occurs from this phase. We 
start with a simulation cell containing only the h.c.p. solid, with the basal plane parallel to one 
face of the cell, and the system is allowed to thermalize at a temperature where coexistence 
is expected. At some instant of time, the simulation is stopped, and the atoms in one half of 
the system are held fixed; the boundary between the two halves is taken parallel to the h.c.p. 
basal plane. The atoms in the other half are raised to a very high temperature, and dynamical 
evolution of these atoms is allowed to proceed so that melting occurs. With this half of the 
system molten, its temperature is now reduced to the original value, the atoms in the other half 
still being fixed. Finally, the atoms in the fixed half are given thermal velocities and released, 
and the whole system is allowed to evolve freely. The system is monitored by calculation of the 
average number density in slices of the cell taken parallel to the boundary between solid and 
liquid. As we shall show, the density in the solid part is a periodic function of slice number, 
while in the liquid it fluctuates rather weakly about its average value. The total energy, 
temperature and pressure are, of course, monitored throughout the simulation. A feature of 
the (JV, V, E) method is that at equilibrium the stress in the solid phase will generally not be 
hydrostatic, and manual adjustment of the cell parameters is needed to achieve hydrostatic 
stress. We shall show in Sec. f| that this does not present a problem. 

2.4 Correcting the coexistence approach 

Both the free-energy and the coexistence approaches are subject to errors, which need to be 
assessed and corrected for. It may be too costly to perform all the ab initio calculations with 
the precision needed to obtain accurate melting properties, so that there are errors due to 
inadequate fc-point sampling, incompleteness of the basis set, or other approximations. In the 
coexistence approach, inevitable differences between the reference and ab initio total energies 
will create errors in the melting properties, which need to be evaluated. The correction of 
errors in the free-energy approach has been extensively discussed elsewhere []|, 17], so our main 
concern here is correction of the coexistence approach. 

Assuming the ab initio energy is calculated with adequate precision, the main error comes 
from differences between U\j and U re {. The key question posed here is therefore: how are 
predicted melting properties changed by small changes in the total energy function? (System- 
size errors will be discussed separately later.) Related questions have been discussed before 
(see e.g. Ref. |14}| ), so we give only a brief summary of the main points. 

The difference Uai — U ie { is denoted by AU, as in the free-energy approach. The change 
of any quantity resulting from the replacement of U Ie f by Uai will also be indicated by A; 
for example, the shift of melting temperature at a given pressure is called AT m . The melting 
temperature is shifted because the liquid and solid Gibbs free energies G l (p, T) and G s (p, T) , 
and hence their difference G ls (p,T) = G l {p,T) - G s (p,T), are shifted. Working at the given 
pressure, we take the variable p as read and express the ab initio value of G ls as: 

G l A l (T) = G l T s cf + (AG ls (T) , (4) 

where the parameter £ is introduced so that the reference melting temperature 7™ can be 
written as a power series: 

T^ = T- f + <X + C 2 ^ + ... . (5) 



Since the Gibbs free energies are equal in the two phases, this T^ 1 is the solution of G l ^\(T) = 0, 
which is: 

G l r s ei (r- f + (X + C 2 ^ . . .) + (AG ls (T- f + cr; + c 2 K + ...)= o . (6) 

Expanding in powers of £ and equating powers, one obtains formulas for T^, T^, etc: 

T' = AG ls (T" ci ) /Sl s 



m — \ J m y / °ref 

m^pjref / 



J m/ J m — ~~FhT J m p,ref/ z - f m T ; {< ) 

^ref 



where S^ f = 5' ef — S^ef is the reference entropy of fusion, C l p s ref = C l p ref — C p ref is the liquid- 
solid difference of the constant-pressure specific heats, and AS ls is the shift of the entropy of 
fusion. We note particularly the implication of the formula for T^. Since entropies of fusion 
are on the order of k-Q per atom, then a difference AG ls of 10 meV/atom implies a shift of 
melting temperature of ca. 100 K, so that substantial errors will need to be corrected for 
unless the reference total energy function matches the ab initio one very precisely. Although 
we have included a formula for T^, it seems unlikely that this will be used in practice, except 
perhaps for a rough estimate of T^/T^, since AS ls would be difficult to compute without 
rather extensive free-energy calculations. 

If simulations are done in the isothermal-isobaric ensemble, then for closely matching Uai 
and U Te {, with small fluctuations of AU, the Gibbs free energy shifts AG 1 and AG S can be 
evaluated using the well-known expansion: 

AG = (AU) Iet - ^(3(5AU\ ci + . . . , (8) 

where 8AU = AU — {AU) Te {, and the averages are taken in the reference ensemble. The ab 
initio simulations presented later were performed in the isothermal-isochoric ensemble, so that 
the quantity that is readily evaluated is AF(V,T), the change of Helmholtz free energy when 
t7 re f is replaced by C/ai are constant volume V. This AF is given by the same formula Eqn (ft), 
but with the averages evaluated in the isothermal-isochoric ensemble. In this case, we need to 
consider the relation between AG and AF, which is readily shown to be: 

AG = AF - Vk t A p 2 , (9) 

where kt is the isothermal compressibility and Ap is the change of pressure when U re { is 
replaced by Uai at constant V and T. 



2.5 Size effects 

Whichever approach is used, coexistence must be treated in the thermodynamic limit. In the 
free-energy approach, almost all the free energy is that of the reference system, for which size 
errors can be made negligible by doing simulations on very large systems. Appreciable size 
errors remain only in the difference -Fai — F re f between the ab initio and reference systems. 
These can only be assessed and corrected for by explicitly calculating AF = AF = Fai— F re f for 
systems of increasing size, as we report elsewhere for the cases of Al || and Fe ||]. Similarly, 
in the coexistence approach, explicit coexistence simulations are performed with the model 
system on very large systems, so that size errors are made negligible. However, the shift of 
melting temperature due to the difference Uai — U Ie f requires calculations of AF of exactly 
the same kind as are needed in the free-energy approach. This means that the size errors are 
essentially the same in the two approaches. We give a quantitative assessment of size errors in 
the melting properties of Fe by the coexistence approach in the next Section. 



3 The melting properties of iron 



3.1 Technical details 

Our coexistence calculations on high-pressure Fe use precisely the same DFT techniques used 
in our free-energy work || [17]], so we give only a brief summary here. The exchange-correlation 
functional E xc is the generalized gradient approximation known as Perdew-Wang 1991 pj, [2q], 



We use the projector-augmented- wave (PAW) implementation of DFT |2^, 27, |2q], a technique 
that shares the properties both of all-electron methods such as full-potential linearized aug- 
mented plane waves (FL APW) |2{J and the ultrasoft pseudopotential method . The calcu- 



lations were done using the VASP code [31, 32]. Details of the core radii, augmentation-charge 



cut-offs, etc. are exactly as in our PAW work on liquid Fe [28]. Our division into valence and 
core states is also the same: the 3p electrons are treated as core states, but their response to 
the high compression is represented by an effective pair potential, with the latter constructed 
using PAW calculations in which the 3p states are explicitly included as valence states. 

The reference model for our coexistence simulations is the embedded-atom model (EAM) 
recently used by Belonoshko et al. || to calculate the high-pressure melting curve of Fe. This 
EAM has the standard form [^] , in which the total energy E^ = J2i &i is a sum °f energies 
Ei of atoms i, with each E$ = E™ p + F(pi) consisting of two parts: first, a purely repulsive 
energy E^ represented as a sum of inverse-power pair potentials E^ = ^(a/rij) n , where 
Tij is the interatomic separation and the sum excludes i = j; second, an 'embedding' part 
E(pi) which accounts for the metallic bonding mainly due to partial filling of the ci-bands. The 
embedding function F(p) is represented as — eCp 1//2 , and the density pi for atom i is given by 
the sum over neighbours pi = J2'j{ a / r ij) m - The parameters in this EAM were determined in 
Ref. H by fitting to full-potential linearized muffin-tin orbital (FPLMTO) energies for typical 
configurations of liquid iron, and they are: n = 8.137, m = 4.788, e = 0.0173 eV, a = 3.4714 A, 
and C = 24.939. 

The aim of our calculations is not to generate the entire melting curve, but to obtain a 
point on this curve at a pressure close to the value of 330 GPa at the boundary between the 
Earth's solid inner core and liquid outer core plj. Our main simulations have been done on 
cells containing 8000 atoms, constructed as a 10 x 10 x 40 h.c.p. supercell. The long axis 
is perpendicular to the basal plane of the crystal, and we refer to it as the z-axis. We have 
tested the adequacy of this system size by performing coexistence calculations on systems of 
up to 20,480 atoms. Within the statistical errors, we were unable to detect any difference 
between results for coexistence pressure and temperature with this system size and those with 
the 8,000-atom system. The mean volume per atom for the calculations reported here was 
v = 7.12 A 3 , which is near the volumes of the solid and liquid in the pressure region of 
interest. Since our calculations are performed at constant volume and constant cell shape (see 
Sec. ^), we checked that non-hydrostatic stresses do not affect the melting properties of the 
system. To do this, we performed the simulations with the two c/a values 1.64 and 1.66, which 
give the P zz diagonal component of the stress tensor slightly larger or smaller than the P xx 
and P yy components in the two cases. We find that the effect on the melting temperature is 
undetectably small. 



In preparing the coexisting solid and liquid as described in Sec. 2.3, the system was initially 
equilibrated at 6100 K, and the high temperature used to ensure complete melting of half the 
system was 5 x 10 4 K. Once melting had been achieved, the liquid part was re-equilibrated at 
6100 K before free evolution of the whole system was started. 



3.2 Results 

In Fig. |l], we show the temperature of the system as a function of simulation time for c/a = 1.64. 
We also show the three diagonal components of the stress tensor; the off-diagonal components 
fluctuate around their average value of zero, so there is no shear stress on the cell. After 



an initial equilibration period, one sees that the temperature and the pressure settle around 
the values T = 6550 ± 100 K and p = 305 ± 1 GPa. In Fig. ^, we display the density profile, 
calculated by dividing the simulation cell into 400 slices parallel to the liquid-solid interface and 
counting the number of atoms present in each slice. The profile shown corresponds to the last 
configuration in the simulation with c/a = 1.64; a similar profile is observed in the simulation 
with c/a = 1.66. It is easy to identify the solid and the liquid regions in the system: in the 
solid region the density is a periodic function of slice number, but in the liquid it fluctuates 
randomly around its average value. This confirms that we do indeed have coexisting solid and 
liquid in the cell. At 305 GPa, the melting temperature reported by Belonoshko with precisely 
the same EAM is 6680 K, which agrees with our value within the combined statistical errors. 

We now correct the EAM reference melting temperature to obtain the fully ab initio melting 
temperature using the techniques presented in Sec. 2.4. In order to assess possible system-size 
errors, the free-energy corrections of Eqs. (|8]J|) were calculated with systems of 64, 150 and 
288 atoms. To do this, we used the EAM reference model to generate long simulations for the 
solid and liquid separately. From these simulations we extract typically 50 and 100 statistically 
independent configurations for solid and liquid respectively, for which we calculate the DFT 
total energies U\i. The differences AU = Uai — U re { are then used to compute the free-energy 
corrections and hence the shift of melting temperature (see Eq. (|7])). The temperature was set 
equal to the value of 6550 K that emerges from the coexistence simulations, and volumes per 
atom of v s = 7.05 and vi = 7.218 A 3 /atom were also deduced from the coexistence simulations; 
the pressure in both cases was 303 GPa. (The very small difference from the pressure of 305 GPa 
mentioned above is of no consequence.) From these two volumes we can extract the volume 
change on melting, which is AV ~ 2.5%, the same value reported as by Belonoshko et al. |J. 
(The value from our ab initio free-energy work was 1.8 % at 305 GPa ||.) We also calculated 
the entropy of melting from the relation TS\ s ei = E\ s ei + pV^ f ; the energy change on melting 
E\ s et was calculated using the two separate simulations of solid and liquid at p = 303 GPa 
and T = 6550K, and we obtained S[ s ei = 0.88 kg /atom; a similar value of S ls can be deduced 
from the work of Belonoshko et al. Q by using their value for V ls and the Clausius-Clapeyron 
relation on the reported melting curve. (The value from our ab initio free-energy work was 
1.07 fee/atom at 305 GPa ||.) The DFT energies were carefully checked for electronic Appoint 
errors, as in our free-energy work H |l7|| . As expected from that work, T-point sampling for 
the 64-atom system underestimates the fully converged energies of liquid and solid by ca. 10 
and 50 meV/atom respectively. With 150 atoms, T-point sampling gives negligible errors for 
the liquid, but an error of ca. —8 meV/atom for the solid. With 288 atoms, we have used 
only T-point sampling, but the indications are that /c-point errors should be negligible for both 
phases. 

We report in Table [j] the results for (A[7) re f and ((<5A[/") 2 ) re f for the three system sizes, 
the results being already corrected for electronic /c-point errors. An important feature of 
these results is that there is no discernible system size effect on these corrections within the 
statistical errors of 5 meV or less, so that systems of 64 atoms are adequate for calculating 
the corrections. From Eq. (^), this imples that the shift of T m due to size errors will not more 
than ~ 50 K. In order to obtain the corrections AG to the Gibbs free energy, we need to 
include the term in Ap 2 in Eq. We find that the pressure differences between the EAM 
and ab-initio systems are only 22 and 12 GPa for solid and liquid respectively, which give this 
term values of ca. 5 meV. From Table [l], we see that the free-energy differences between the 
ab initio and EAM systems have the effect of stabilizing the liquid with respect to the solid by 
about 35 meV/atom compared with the EAM. In order to obtain the correction to the melting 
temperature from Eq. (|7j), we use the value for the entropy of melting S l r s ei of 0.88/cB/atom 
quoted above. This gives the first-order correction = —450 K, so that our corrected T m at 
305 GPa is 6100 ± 100 K. This result is in very close agreement with the free-energy approach, 
which at p = 305 GPa gives T m = 6100 K. (This value is somewhat lower than the preliminary 
result of our free-energy calculations jl ; the downward revision came from a careful reanalysis 



of the anharmonic free energy of the h.c.p. crystal, as reported in Refs. [||, |l7|.) 

Although the EAM of Belonoshko et al. |8| mimics our ab initio systems reasonably well, 
we have found that the model can be still further improved by refitting it to our ab initio 
energies of the solid and liquid. In doing this, we allowed only the repulsive potential and 
the strength of the embedding energy to change, so that only the parameters n, e and G are 
allowed to vary. These were adjusted to minimise the 5AU fluctuations for the solid, while also 
maintaining the correct pressure. The new parameters are n = 5.93, e = 0.1662 eV, C = 16.55, 
the parameters a = 3.4714 A and m = 4.788 retaining their original values. We then repeated 
the coexistence simulation using this new EAM and obtained a coexistence temperature of 
T m = 6200 ± 100 K at a pressure of p = 323.5 ± 1 GPa. The very small free-energy corrections 
(Table slightly stabilize the solid, resulting in an increase of melting temperature of ca. 50 K 
to give 6250±100 K. The differences with the ab-initio pressures are negligible in this case, and 
the term in Ap 2 in Eq. @ does not contribute to the free energy difference. At this pressure, 
the melting temperature from our free energy approach was T m = 6290 K [||, so that once 
again we find very good agreement between the two approaches. 

4 Discussion and conclusions 

We have advocated the aim of calculating the melting properties that follow from a chosen 
approximation for the exchange-correlation functional E xc , with the errors due to all other 
approximations being made negligible. We have stressed that in the free-energy approach the 
total energy of the reference model does not have to agree exactly with the ab initio total 
energy, since this approach includes the calculation of the free energy difference between the 
two systems. Up to now, no allowance has been made for such differences in the coexistence 
approach, but we have shown how corrections can be made so that the final results are again 
independent of the reference model. Our practical results for high-pressure Fe demonstrate 
that the two approaches then give the same melting temperatures, as would be expected. 

Our analysis has implications for the way in which reference systems are constructed. In the 
free-energy approach, the crucial requirement of a reference model is that the fluctuations of 
the difference AU between the ab initio and reference energies be as small as possible; we have 
the freedom to use different reference models for the solid and liquid phases, and a constant 
offset between the ab initio and reference energies is of no consequence. In the coexistence 
approach, it seems clear that the model should be chosen so that the AG corrections, and the 
resulting shift of melting temperature away from that of the reference system, be as small as 
possible. This demands more of the reference system than in the free-energy approach. It is 
necessary but not sufficient that the AU fluctuations be small, in order that AG be small (see 
Eq. @). In addition, a constant energy offset will shift T m , if the offset is not identical in the 
two phases, since the same model must simultaneously reproduce the energetics of both. This 
implies that force matching |?], |l9| may not always be a reliable way of constructing reference 
models: even if the reference and ab initio forces match precisely in both phases, there may 
still be a non-cancelling AU offset that will shift T m . In this sense, energy matching is safer 
than force matching. 

Concerning finite-size errors, we have shown that these affect appreciably only the free- 
energy differences between the ab initio and reference systems, since in both approaches the 
calculations on the reference model will always be done on systems large enough to make 
size errors negligible. If anything, these errors may be more troublesome in the coexistence 
approach, because of the need to make the reference model fit both phases at once, so that 
the free-energy corrections may sometimes be larger. However, in our coexistence calculations 
on high-pressure Fe, we have shown that these errors are unlikely to shift T m by more than 
~ 50 K. 

In commenting on the strengths and weaknesses of the two approaches, we note that they 
differ mainly in the way of treating the thermodynamic properties of the reference model. In 



both approaches, corrections are then needed to obtain the properties of the ab initio system, 
and these involve ab initio simulations to perform thermodynamic integration or to compute the 
free-energy differences of Eq. (@). The melting properties of the reference model are probably 
more straightforward to calculate by the coexistence approach, since the free-energy approach 
requires rather intricate thermodynamic integrations. However, the heaviest computational 
effort comes in the calculation of the free-energy corrections, and here we believe the free- 
energy approach may have the advantage, since the effort can be reduced by using different 
reference systems for the two phases. In practice, we have found it very helpful to use both 
approaches, as in the work on Fe reported here. 
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(AU) 


ref/iV 


{{5AU) 2 


)rcf/iV 


N 


Liquid 


Solid 


Liquid 


Solid 


64 


-6.940 ± 0.003 


-6.909 ± 0.002 


0.023 ± 0.004 


0.014 ± 0.003 


150 


-6.934 ±0.001 


-6.912 ±0.001 


0.023 ± 0.003 


0.009 ± 0.002 


288 


-6.939 ± 0.001 


-6.909 ± 0.001 


0.024 ± 0.004 


0.014 ± 0.003 



Table 1: Thermal average (Af7) re f of the difference AU = Uai — U rc { of ab initio and reference 
energies, and thermal average ((<5A£/) 2 ) re f of the squared fluctuations of 5AU = AU — (Af7) re f, 
with averages evaluated in the ensemble of the reference system and normalized by dividing by 
the number of atoms N (eV units). Results from simulations with difference N are reported 
for liquid and solid Fe for the thermodynamic state p = 303 GPa, T = 6550 K. Reference 
model is the embedded-atom model of Belonoshko et al. ||]. 



N 



(i I 



(AU) IC{ /N ((5AU) 2 ) TC{ /N 
Liquid Solid Liquid Solid 



7.213 ±0.002 7.200 ±0.002 0.018 ± 0.003 0.009 ± 0.002 



Table 2: Thermal averages (A[7) re f and ((5AU) 2 ) m f (eV units, see caption of Table 1) for the 
EAM reference model obtained by refitting that of Belonoshko et al. J8|] . 
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Figure 1: Time variation of temperature (upper panel) and the three components of stress 
tensor P xx (solid curve), P yy (dashed curve) and P zz (dotted curve) (upper panel), during 
a simulation of solid and liquid Fe coexisting at a pressure of 305 GPa. Simulations were 
performed on a system of 8000 atoms using the embedded-atom potential of Belonoshko et 
al. |§, with c/a ratio of the h.c.p. solid equal to 1.64. 
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Figure 2: Density profile in a simulation of solid and liquid Fe coexisting at a pressure of 
305 GPa. The system is divided into slices of equal thickness (0.35 A) parallel to the solid- 
liquid interface, and graph shows number of atoms in each slice. Simulations were performed 
on a system of 8000 atoms using the embedded-atom potential of Belonoshko et al. ||. 



